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ABSTRACT 


Numerical methods for solving the flow equations in cylindrical or spherical coordinates 
should be able to capture the behavior of the exact solution near the regions where the 
particular form of the governing equations is singular. In this work we focus on the 
treatment of these numerical singularities for finite-differences methods by reinterpreting the 
regularity conditions developed in the context of pseudo-spectral methods. A generally 
applicable numerical method for treating the singularities present at the polar axis, when non- 
axisymmetric flows are solved in cylindrical coordinates using highly accurate finite- 
differences schemes (e.g., Pade schemes) on non-staggered grids, is presented. Governing 
equations for the flow at the polar axis are derived using series expansions near r=0. The 
only information needed to calculate the coefficients in these equations are the values of the 
flow variables and their radial derivatives at the previous iteration (or time) level. These 
derivatives, which are multi-valued at the polar axis, are calculated without dropping the 

accuracy of the numerical method using a mapping of the flow domain from (0,R)*(0,2rc) to 

(-R,R)*(0,7t), where R is the radius of the computational domain. This allows the radial 

derivatives to be evaluated using high-order differencing schemes (e.g., compact schemes) at 
points located on the polar axis. The proposed technique is illustrated by results from 
simulations of laminar-forced jets and turbulent compressible jets using large eddy simulation 
(LES) methods. In term of the general robustness of the numerical method and smoothness 
of the solution close to the polar axis, the present results compare very favorably to similar 
calculations in which the equations are solved in Cartesian coordinates at the polar axis, or in 
which the singularity is removed by employing a staggered mesh in the radial direction 
without a mesh point at r=0, following the method proposed recently by Mohseni and 
Colonius [1], Extension of the method described here for incompressible flows or for any 
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other set of equations that are solved on a non-staggered mesh in cylindrical or spherical 
coordinates with finite-differences schemes of various level of accuracy is immediate. 

Key words: coordinate system singularity, polar axis, turbulence, large eddy simulation 
method, cylindrical system, collocated grids, finite differences method 
Numerical analysis category: 65P10, 65P99 
Genera] category: 76Aaxx 
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1. BACKGROUND 


The singularities at the centerline of a cylindrical coordinate system are due to the 
presence of terms containing the factor 1/r, where r is the radial distance, in the equations 
governing the flow. The flow field itself does not have any singularity at the polar axis. For 
axisymmetric flows symmetry conditions may be used to remove these singularities, but in 
the simulation of non-axi symmetric flows this problem must be faced. Furthermore, as the 

computational domain is defined as (0,27t)*(0,R) one has to specify numerical boundary 

conditions at n=0, even if physically there is no boundary at the polar axis. 

Several numerical methods have been proposed to address the singularity of the flow 
equations in cylindrical or spherical coordinates. At first sight they seem to vary greatly 
depending on whether a pseudo-spectral, finite-volume or finite-differences framework is 
adopted, but in fact there are many similarities among them. This is especially true for the 
new finite-difference method on non-staggered grids proposed here and the treatment used 
with pseudo-spectral methods. There are also common elements with methods that use 
l’Hopital’s rule [2] or use a ‘shifted’ distribution of points in the radial direction [1] to 
eliminate the points on the polar axis. 

The main idea behind using spectra] methods in cylindrical coordinates is to seek an 
approximation using polynomial expansions in the radial direction that satisfy some regularity 
conditions so as to insure a well-behaved solution near the polar axis. O’Sullivan and Breuer 
[3] studied the growth of linear and non-linear disturbances in pipe flow. They solved the 
Navier-Stokes equations using Fourier transforms in both the streamwise and azimuthal 
directions, and Chebyshev polynomials in the radial direction. The singularity at the polar 
axis was avoided by using information in wavenumber space related to the form of the series 
expansions of the velocity and pressure variables at r=0, and by mapping in the radial 
direction to increase resolution near the polar axis. Patera and Orszag [4] in a related study 
used a similar spectral method but they employed parity relations for the Chebyshev 
coefficients in the expansions to remove the coordinate singularity of the cylindrical system. 
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Fully developed pipe flow DNS calculations in cylindrical coordinates using spectral methods 
were reported in [5] and [6]. The method in [5] used Jacobi polynomials as the expansion 
basis in the radial direction close to the polar axis and Legendre-Lagrangian interpolants away 
from the axis, while the latter used B-splines polynomials expansions. An adaptation of the 
method proposed in [5] for spherical coordinates, where one is faced with essentially the 

same problem at the two singularity axes corresponding to polar angles <p=0 and (p= 7 t can be 

found in [7]. Zhang et al. [5] were able to recast the governing system of equations in the 
wavenumber space so that the singularities due to coefficients of the form 1/r and 1/r 2 were 
removed. This was essentially done by taking advantage of the particular form of the series 
expansions in wavenumber space of the pressure and velocity components. However, they 
still had to use a special form of Jacobi polynomials in the radial direction close to the origin 
in conjunction with 1 Hopital’s rule to remove potential singularities arising because of terms 
where both the numerator and the denominator tend to zero at approximately the same rate 
near the polar axis. Loulou [6] followed a similar approach, but the regularity conditions in 
terms of having a certain behavior in r near the origin in the polynomial expansions of the 
Fourier modes were imposed as constraints on the coefficients of the B-splines. These 
regularity conditions are essentially the same ones we are going to use in our method, but in 

the present work everything is formulated in the physical (r,0,z) space as opposed to 

wavenumber space (r.lq,,!^). In the spectral methods the essential information is the order of 
the leading term in the polynomial expansions in r for each mode of the Fourier expansion in 
0. No explicit equations at the polar axis are derived, and there is no need to calculate the 

values of the coefficients in these expansions. Finally, a somewhat different approach was 
proposed in [8] in the context of pseudo-spectral (interpolatory) methods. Using the 
governing differential equations pole conditions were derived in [8] which were subsequently 
used as numerical boundary conditions at the coordinate singularity. The advantage of the 
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method is that these pole conditions are developed in physical space and are easier to 
implement into a standard pseudo-spectral method. The method was illustrated for solving a 
Poisson equation on the unit disk. 

Finite- volume and finite-differences methods are less accurate than spectral methods, but 
they can handle complex geometrical configurations more easily. Finite-volume methods for 
flows in cylindrical geometries were employed by Eggels et al. [9], Akselvoll and Moin 
[10,1 1] and Boersma et al. [12], among others. The main difficulty here is related to the fact 
that the azimuthal and streamwise fluxes are not defined at the centerline. Due to the 
definition of velocity components at the center of the volume cells, the singularity for the 
azimuthal and streamwise fluxes is removed. However, if the control volume closest to the 
centerline is not wedge shaped, the flux at the first control-volume surface off the polar axis 
whose normal is oriented in the radial direction must be evaluated. In one of the 
implementations of their numerical method, Eggels et al. [9] used first- and second-order 
one-sided differences to extrapolate the radial velocity on the closest faces to the centerline. 
A similar procedure was used to calculate the diffusive flux at these locations. Loss of 
accuracy near r=0 was found not to be very significant due to the relative smooth DNS 
velocity fields in that region for the pipe flow. They also tried to use wedge-shaped elements 
at the centerline, the advantage being that the singularity is automatically removed as the 
radial flux is identically equal to zero on the faces ‘located’ at the centerline. In this 
implementation the method was found to give unrealistically high values for the r.m.s. 
fluctuations near the centerline. However, one should point out that Boersma et al. [12], 
using a similar method, did not observed such problems in their DNS simulations of 
incompressible jets. Akselvoll and Moin [10] used a finite- volume method to solve for the 
turbulent flow in a coaxial jet combustor using LES. They used regular-shaped elements 

near the centerline, with the closest faces located at Ar/2 from the singularity axis. Provided 
that the radial and azimuthal velocity components u r and u e are well defined at the centerline 
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for all 0’s (they are multi-valued at r=0), the calculation of the fluxes on these faces presents 
no special problem. They interpolated the values for the radial velocity at every azimuthal 
location 0 on the centerline using the corresponding values at Ar/2 for two azimuthal angles, 

0 and Q+n, and accounting for the change in sign for the radial velocity across the centerline 

in the limit r-»0. A similar treatment was followed for the azimuthal velocity at the 
centerline that was needed to compute some viscous fluxes in the momentum equation for the 
azimuthal fluxes. 

Verzicco and Orlandi [13] developed a second-order finite-differences scheme in 
cylindrical coordinates. The main feature of their method was the introduction of a radial flux 
(ru r ) on a staggered grid. As for finite- volume methods, this is the only flux to be calculated 
right at r=0, and its value there is obviously equal to zero. The DNS study in [14] used this 
algorithm to calculate the fully turbulent flow in pipes with rotating walls. Their method 
seems to be very appealing in the context of second-order schemes, but extension to higher- 
order schemes is not straightforward. 

Griffin et al. [2] used finite differences and L’Hopital’s rule to recast the governing 
equations at the origin. For all the terms involving coordinates singularities they calculated 
the radial derivatives using one-sided second-order-accurate finite differences, while the 
azimuthal derivatives for the variables that are multi-valued at the origin were formally 

evaluated by using the values at neighboring azimuthal locations and A0=27t/N e , where N e is 

the number of cells in the azimuthal direction. This resulted in a set of new equations at the 
origin that do not contain any singularities. Formally this treatment is not very rigorous as all 
these values are physically located at the same point, but the second-order method was found 
to produce smooth solutions near the polar axis at least in the framework of the inviscid 
calculation of the I.C. engine that they considered. However, its use for DNS or LES 
simulations using higher-order schemes is doubtful, and in the next section we will discuss 
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in a more rigorous way using series expansions the errors introduced by the use of 
L’Hopital’s rule to remove the singularities at the polar axis. This will also answer a 
question which they raised in their paper, whether or not a fixed limit exists for terms like 

(l/r)dp/30, which they recasted using l’Hopital’s rule intod 2 p/dr00. We show that in this 
case such a limit exits and can be written using Fourier series expansions of the pressure in 
terms of the m=l and m=-l modes (the pressure cross-derivative is multi-valued at the 

origin). For this particular example the Fourier expansions of l/rdp/ 00 and 0 2 p/0r00 give 
the same limit, however we also show other examples (e.g., Laplacian operator that is 
present in the viscous term) where this kind of approximations using 1’Hopital’s rule will 
lead to a wrong result. 

Finally, Mitchell et al. [15], Freund [16] and Boersma and Lele [17] used compact finite 
differences on non-staggered meshes in cylindrical coordinates to compute the radiated jet 
noise of compressible jets. In these studies, the equations at the polar axis were solved using 
a Cartesian coordinate system. The directions of the Cartesian system were arbitrarily 

chosen, and the multi-valued variables (u r , u e ) for the other directions were obtained by 

rotation of the system at r=0 where only modes m=l and m=-l may exist. As the present 
numerical scheme is very similar to the ones used in these studies, more details can be found 
in the validation section. Another option to remove the singularities at the polar axis is to 

eliminate the points at r=0 by distributing the points in the radial direction starting with Ar/2 

(for a uniform grid spacing in r) and mapping the domain (0,2 ti)*( 0,R) into (0,7t)*(-R,R) 

when radial derivatives must be computed, such as there is no numerical boundary conditions 
to be specified at the polar axis. The radial derivative stencils will span the centerline without 
evaluation at r=0. This is essentially the method discussed in [1]. We will use this idea but 
recast it for the case in which we have points at the polar axis. As we will derive a set of 
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exact equations which are well defined on the polar axis we will not have to avoid the 
presence of such points at the centerline. 


2. TREATMENT OF THE EQUATIONS AT THE POLAR AXIS 


Let’s suppose that the system of governing equations may be written as: 

^ = RHS(Q) (1) 

dt 

where, for instance, in the case of compressible three-dimensional flows 
Q=(pu x ,pu r ,pu e ,p,e), p is the density, e is the energy and the right hand side term (RHS) 

contains the usual operators in cylindrical coordinates associated with the continuity or the 
transport equations of momentum and energy, including the terms associated with the sub- 
grid scale contributions in the case of LES simulations. 

The numerical method is general enough so as the particular form of the operators in the 
RHS is not really important. The same is true if more equations are included in Eq. (1), such 
as transport equations for passive scalars, or for turbulence quantities in the case of RANS 
calculations. The only important difference between the different equations in Eq. (1) is 
determined by whether the variable in the left hand side of Eq. (1) is single valued or multi- 
valued at the polar axis. Following Boyd [19], the most general expansion of any function 
that is not multi-valued at the polar axis (scalars, pressure, streamwise velocity) can be 
written as follows: 


F(r, 0) = I (f m (r) • cos(m0) + g m (r) • sin(m0)) 

m-0 


( 2 ). 
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where f m and g m are polynomials in r that have m-th order zeroes at r=0. If m is even, f m and 
g m are both symmetric around r=0 and their power series contain only even powers of r, 
while if m is odd, f m and g m are both antisymmetric and accordingly their power series 

contain only odd powers of r. A similar statement holds for ru r and ru 9 which takes care of 

the expansion form for multi-valued variables at the origin. These conditions imply that the 
most genera] series expansion of a single valued quantity (S) at the polar axis can be written 
as: 


S(r,0)= I r 

m=0 


m 


( 

ZCt : 

vn=0 


2n 


ran 1 


■cos(m0) + X r 

m=0 


m 




_2n 


• sin(m0) 


\n=0 


(3) 


while the expressions for multi-valued quantities (e.g., u r and u 9 ) assume the following form: 


1 


M(r,0) = -£A On r 2n + £ r 
r n=l m=l 


m-1 


ZA^r 2 " I • cos(m0) + £ r 


m-1 


\n=0 


2B m „r 2 " 


m=l 


\n=0 


■ sin(m0) 


J 


(4 > 

A sketch of the proof of Eq. (3) can be given starting with the most general Fourier 
expansion of S in the form of Eq. (2) and requiring that all the terms are regular at the origin. 
The Fourier expansion in terms of complex variables can be written as: 


F(r,0)= l(f m (r)e i "‘ e ) ( 5 ) 

m=0 


where f m (r) is a polynomial in r with complex coefficients. Each term in the summation can 
be written as: 
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( 6 ) 


f m (r)e"" e = ic n (m)rV m9 = r”V " 9 |c„(m)r”- m = w m f c„(m)r ] 

n=0 n=0 n=0 

where w=y+iz=re' 0 . As w m is regular at origin for all m>0, the requirement that each term 
in Eq. (5) is regular implies c n (m)=0 for n<m, as 1/r 5 is not regular at origin for s>0. Also 

r s = (V(y 2 +z 2 )) is regular only for even values of s. For odd values of s the derivative of 

r 5 with respect to y or z will contain 1 / Vr which is not analytic at origin. 

The form of the series expansions for multi-valued quantities (e.g., for u,) can be 
deduced by observing that u r can be written in terms of the Cartesian velocity components u y 
and u z as: 


u r =u y cos(0)+u z sin(0) 


(7) 


where the form of u y and u z ’s expansions is given by (3). By regrouping the terms and using 
algebraic identities of the form 2cos(0)cos(m0)=cos((m-l)0)+ cos((m+l)0) one can easily 
recover Eq. (4). 

As we previously mentioned any scalar or Cartesian velocity component is uniquely 
defined at the origin, so one can write: 



= 0 


( 8 ) 


This relation holds in particular for u z — u r sin(0)+u e cos(0), so by taking the derivatives 


with respect to 0 one can write: 
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( 9 ) 


0 = 


9u r 

de~ 


U 0 




cos(0) 


r=0 


As the Eq. (9) should hold for any 0, one obtains: 


9u r , 9 u b 

— = u e and — = -«, a.r=0 


( 10 ) 


There is another important constraint on the general form of the series expansions for u r and 
u 0 • If A«Bjp,AW and are the coefficients of the series expansions for u r and u e in 

Eq. (4), the following relation holds for all i>l: 


Aff-Bff a„dBg> = -A!;> 


( 11 ) 


This relation can be obtained by using the corresponding series expansions given by Eq. (3) 
for u y and by Eq. (4) for u r and u e , respectively. These series expansions are plugged into: 


u y =u r cos(0)-u 9 sin(0) (12) 

and using algebraic identities similar to the ones employed to deduce the general form of the 
series expansions for multi-valued variables at the polar axis (Eq. 4), one can see that the 

polynomial multiplying cos(20) in the RHS of Eq. (12) has the form: 


12 


^(bJo* + Ajo + r 2 (—) + — ) 


(13) 


As this polynomial should have a double zero at the origin, one should require B,^ = -A(«\ 
The other relations in Eq. (11) follow similarly from analysis of the form of polynomials 
multiplying the other modes in the RHS of Eq. (12). 

By calculating the derivatives with respect to 0 and r of the series expansions given by 

Eqs. (3) and (4) for all operators present in the RHS of the governing equations Eq. (1) and 
taking the limit r-> 0 a new form of the governing equations that is valid at r=0 is obtained. 
These are a set of exact equations at the polar axis, provided that we can calculate exactly the 

coefficients A^, B mn , a ma , for all terms present in the RHS of Eq. (1). However, for a 

system of PDE’s with second-order radial derivatives, as is the case for the Navier-Stokes 
equations, it is sufficient to calculate at most the coefficients whose indices m and n vary 
between 0 and 2, as we show next. The final expressions in the RHS of Eq. (1) are 
dependent on the particular set of equations that is solved (inviscid, laminar, turbulent, 
compressible, incompressible, etc.). However, the present method is valid for any system of 
equations that can be written in cylindrical (or spherical) coordinates in the form of Eq. (1). 
Hence, reference is made to the Appendix for the final expression of the RHS corresponding 
to the compressible turbulent flow equations that were used in the present flow calculations to 
validate this method. In the following discussion we will focus our attention on pointing out 
some general features. 

At this point it is relevant to look at the expressions for the series expansions of a single- 
valued variable at the origin, let’s say the streamwise velocity u x (coefficients in series 

expansions are labeled and (3^), and for a multi-valued one, the radial velocity u r 

(coefficients are A^ andB^), as well as their first and second radial derivatives, in the 
limit r— >0: 
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(14) 


u x = «oo 


= aSo cos(0) + Pi? sin ( 0 ) 

3 2 u 


dr 2 


h = 2a?, ) + 2a ( 2 ? cos(20) + 2p$ sin(20) 


u r = A$ cos(0) + B$ sin(0) 

^ = A{? + A ( 2 g cos(20) + B$ sin(20) 
or 


3 2 u r _ 


dr 2 


= 2Aj r 1 ) cos(0) + 2BS r , ) sin(0) + 2A$ cos(30) + 2Bg sin(30) 


(15) 

(16) 

(17) 

(18) 

(19) 


As expected the series expansions of scalar quantities (e.g., u x ) contain only the m=0 mode, 
while those of u r and u 9 contain only the m=l (cos(0)) and m=-l (sin(0)) modes. 

Meanwhile, the first derivative of a scalar quantity is multivalued as it contains the m=l and 
-1 modes, while the second derivative contains the m=0, 2 and -2 modes. In particular, this 
applies to the pressure, p, which is a single-valued scalar, but dp/ 8r is multi-valued at r=0. 
If we go back to Eq. (1) for the u x variable, and take the limit as r — > 0 , we obtain: 


0u x = faro 

dt at 


= RHS 


( 20 ) 


where RHS U> contains terms coming from the convective and viscous operators that would 

obviously involve modes m=0,l,-l,2,-2,3 and -3 (see Appendix). In order for Eq. (20) to 
be valid, the coefficients resulting from the series expansions of the different variables and 
their derivatives should be such that the coefficients multiplying all but the m=0 mode in 
RHS U must be equal to zero. Indeed this is what happens when we sum the contributions 
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coining from all the terms contained in RHS U and use Eq. (11). However, individual terms 

may contain non-zero coefficients multiplying the m^O modes. In fact, in can be proven that 
this property holds not only for the RHS of the equations describing the transport of a single- 
valued variable at the polar axis as a whole, but also for independent operators such as 
convective, viscous, dilatation, viscous dissipation in the energy equation, etc. For instance, 
the dilatation operator, which is a scalar quantity and thus should contain only the m=0 
mode, can be expressed as: 


d u x , for 

dx dr 


/ du f 
u r + • 


v 


ae 


^«L + 2A ( '> 
dx +2A °' 


( 21 ) 


where the streamwise derivative can be calculated with exactly the same method used (6 th 
order compact differences in the present work) for points situated away from the polar axis. 
Same observation is true for the continuity equation, recasted as an evolution equation for the 
density: 


dp _ f dpu x | 1 dipu r 1 apue ^j _ | 

dt dx r dr r dd J 


f da (p) rr (x) 

g «00 a 00 + A (r) oc (p) 4- R( r >ft(P> + 7 

+ A 10 a 10 +tJ ioPio +ZA Ol a OO 




dx 


J 


(22) 


In Eq. (22) the series coefficients for the density p are andp^. Another example is the 


Laplacian operator (e.g., the viscous term in the u x -momentum equation, or more exactly the 
term corresponding to the constant viscosity part): 


a 2 u x a 2 u x i du x i a 2 u x _ 

dx 2 dr 2 r dr r 2 d0 2 


d 2 a 


dx 2 


(x) 

00 4a[ ) x 1 ) 


(23) 
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For the u r and u 0 -momentum equations the RHS will contain only the m=-l and m=l 

modes, once contributions from all the terms are added. These two equations are identical in 
the limit r — > 0, thus it is sufficient to calculate the limit for only one of the equations, e.g.: 


du r _ 0Aj? cos(0) + B$ sin(0) 

at " at 


= RHS m_1 cos(0) + RHS m=_1 sin(0) 


(24) 

In fact, we will end up with two scalar equations corresponding to the m=l and -1 modes, 
respectively as one can see from the previous relation, the reason being that as r -» 0, u r can 

be obtained from u e by a counterclockwise rotation of 7i/2. This means that if the radial 

component at the polar axis is given by Eq. (17), the azimuthal component can be written as: 


u e = A So cos ( 0 ) + B io } sin ( 0 ) = -AS? sin(0) + B$ cos(0) (25) 


where use was made of relation (11). Relation (23) can also be used to show an inherent 
problem with methods in which l’Hopital’s rule is applied numerically to remove singularities 
at the polar axis (e.g., [2]). In these methods the Laplacian will be calculated at the polar axis 
by numerically evaluating the following expression: 





a 2 

dr 2 




d9 2 


a^og> 

dx 2 


+ 4a? 1 ) - 4 a 2 ? cos(0) - 40$ sin(0) 


(26) 

which introduces a finite error, as the coefficients of the m=2 and m=-2 modes are now 
present. 

The last step needed to complete the presentation of the present method is to describe 
how the asymptotic series coefficients that are needed to evaluate the RHS in Eq. (1) are 
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computed. We will show that all that is required to calculate these coefficients accurately is to 
be able to estimate numerically the first and second order radial derivatives of all the variables 
in RHS with the same order of accuracy as at points away from the polar axis. In particular, 
in our compressible flow solver we are using 6 th order Pade schemes to calculate the radial 
derivatives. To do this, the following algorithm, similar to the one used in [1] is adopted. 
The computational domain is mapped at every x=constant, such as there is no need to specify 

numerical boundary conditions at r=0. The mapping function (r,0)— » (f,0) is: 



0 < 0 < 7t 

0<r<R 


and 


r = -r, 

a for 

0 = 0-71 


71 < 0 < 271 
0<r<R 


(27) 


All variables are identical in both systems of coordinates for O<0<7t. The rules concerning 

the sign changes for the different variables and operators in the subdomain defined by 0<r<R 
and 7i<0<27i follow from Eq. (27): 


u f = 


df 


dr 


’ 


dx 


dx 


?d0 


u- = 

0 dtx 


rd0 


dtx 

a a a 


= -u e , 


a? ar’ a§ a© 


and s = s for any scalar variable s 


(28) 


As all the signs of scalar quantities and streamwise velocities are left unchanged by the 
mapping, Eq. (28) is sufficient to determine the sign of all the terms involving radial 
derivatives in the RHS of Eq. (1) in the mapped domain. The radial derivatives are now 
taken from -R to R with r=0 being a regular interior point instead of a ‘numerical’ boundary 
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point. For example, before taking the radial derivative of (rpu e u r ), which appears in u e 
momentum equation, one has to multiply the values for 0=(7i,27t) by -1, because of the three 

changes of sign (for r, u e and u r ) required by the mapping. One should point out that even if 

a scalar quantity is single valued at origin, its radial derivative is not. Once the radial 
derivatives are calculated, the inverse transformation of Eq. (28) is taken such as the 
azimuthal and streamwise derivatives are estimated with the usual coordinate definition 
(0,R)*(0,27t). 

Next, we will expand on the calculation of the required coefficients in the asymptotic 
expansions for u r , and its radial derivatives. All other variables are treated in a similar 

fashion using the appropriate series expansions. Suppose that N e +1 is the number of points 
in the azimuthal direction (with modes 0,±1, ±2,..., ± N/2). According to Eq. (17), only 
the m=l and -1 modes should be present in a series formed by these N 0 values. We will 
formally associate the N 0 values with 0 going from 0 (N e =1) to 2n (N 0 +l). A Fast Fourier 

Transform can, in principle, be applied and the coefficients of the m= ± 1 modes will be a\q 

and B|q> , the other coefficients should be very small. This can also be done for the first and 
second radial derivatives using their asymptotic expansions in Eqs. (18) and (19) to find 

A 0i’ A 20’ B 20 and A^B^aS.bS, respectively. Once these coefficients are calculated 
the limits of terms involving azimuthal and cross derivatives are also known using Eq. (4) 
and (11), while streamwise derivatives of these coefficients are evaluated with the same 
scheme used for interior points. However, there is a more economical way to calculate these 

coefficients. Supposing that N 0 is divisible by 8, meaning that for eveiy element in the series 
there is another situated exactly n /4 away, one can take advantage of the properties of sin(0) 
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and cos(0) functions to evaluate the coefficients in Eqs. (17), (18) and (19). For instance, 
using Eq. (17) for 0 and 0 +ti/ 2 a system of two linear equations with two unknowns ( Aj^ 
and Bj^) is obtained: 

u r (0) = Ajq cos(0) + Bjq sin(0) 

(29) 

u r (0 + 7t/2) = -A^ sin(0) + Bjq cos(0) 

To eliminate the bias toward a certain direction, one can solve the above system for every 
0=(2k I N 0 ) (n-1) with n=l to N e and average the results to get final values for A$ and B$. 

Same kind of treatment can be applied to determinate A^.A^.B^ by first getting the 

coefficient A^ and then using Eq. (18) for 0 and 0+7t/4 to calculate the remaining two 

coefficients. The coefficients involved in the expression for the second derivative of u r can 
be determined more easily if we try to obtain at the same time the coefficients appearing in the 

expression for the second derivative of u 9 : 
d^u 

' = 2 A i ? cos ( 0 ) + 2B n ) sin ( 0 ) + 2B 30 cos(30) - 2 A^J sin(30) (30) 

where use was made of Eq. (11). By writing Eqs. (19) and (30) at 0,0+7i/2,0+7t/4,0+37t/4 

and solving three systems of two linear equations, one can determinate the values of all 
coefficients similarly to the procedure employed above. 
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3. VALIDATION OF PROPOSED NUMERICAL ALGORITHM 

The development of this new technique was driven by the interest of the authors in 
calculating the turbulence dynamics and acoustic radiation of turbulent compressible round 
jets using LES. This problem is especially relevant to the design of more effective jet 
engines with reduced noise emissions. The quality of the jet-noise predictions is determined 
by the accuracy (resolution) of the numerical method, and by the boundary condition 
(physical as well as numerical at the polar axis) treatment and the quality of the mesh. The 
general numerical method is described in details in [18], while the details of the 
implementation of the dynamic LES model in the original DNS code can be found in [17]. 

The numerical method employs compact six-order Pade schemes [20] for the spatial 
derivatives in the non-homogeneous directions, and Fourier spectra] methods in the 
homogeneous (azimuthal) direction. The number of modes is dropped near the polar axis so 
that the CFL constraint will be determined by the radial (or axial) spacing. No clustering of 
the points near r=0 is needed. The solution is advanced in time using a four-step Runge- 
Kutta method. Zonal boundary conditions with artificial damping are employed near the inlet 
and outlet to absorb outgoing disturbances and to avoid spurious noise generation at these 
boundaries. Non-reflecting boundary conditions are used at the lateral boundary r=R 
(=11R 0 ), as well as damping terms in a layer close to this boundary. The Reynolds number 

is Re=U(2R o )/v=36,000 and the Mach number Ma=U/c„=0.9, where R 0 is the initial radius 
of the jet in the inlet section. In all the calculations discussed here the computational grid 
consists of 192*128*64 points in the (x,r,0) directions, which is about an order of 
magnitude coarser than the grid used by Freund [16] to calculate a similar jet at Re=3,600 
using DNS. The time step is At=O.01R o /U, corresponding approximately to CFL=1. The 

mean flow distribution at the inlet plane is assumed to be a rounded top-hat profile with 
periodic disturbances in the streamwise direction. 
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In a first test case no randomized forcing is applied and the LES model is turned off. A 
sinusoidal perturbation at a Strouhal number St=0.5, corresponding to the most amplified 
wave, is imposed on the streamwise velocity profile. The forced-jet solution is quasi- 
axisymmetrical and laminar as seen from the total vorticity contours in Fig. la. The effects of 
the different treatments at the polar axis are investigated in Fig. 2, where snapshots of the 
dilatation fields in the region x>20R o are shown for different treatments at the polar axis. We 
choose the dilatation because this quantity is very sensitive to the centerline treatment and 
also provides relevant information for the jet acoustics in the near field. Case LI 
corresponds to the method where the equations are solved in Cartesian coordinates, case L2 
corresponds to a grid with no points placed at the polar axis (method described in [1]), while 
case L3 corresponds to the treatment using series expansions at the centerline. In the second 
test case randomized azimuthal forcing is applied in the input plane to trigger the three- 
dimensional instabilities for the LES calculation. All other parameters of the simulations, as 
well as the computational mesh are kept the same as those used in the first test case. As seen 
in Fig. lb the jet undergoes transition to turbulence at the end of the potential core situated 
around x=15R 0 . Results are shown in Fig. 3 for two simulations, the first (Tl) uses the 
method of Mohseni and Colonius [1], while in the second (T2) the series expansion 
treatment is employed. 

Considerable success was reported by Freund [16] for DNS simulations of jets al lower 
Reynolds numbers using the present numerical method with the equations at the centerline 
solved in Cartesian coordinates. As we are interested in simulating jets at very high 
Reynolds numbers (2 to 3 orders of magnitude higher than the ones at which DNS is 
presently possible) the robustness of the method on coarser meshes is an obvious 
requirement. Our simulation showed that even for a laminar forced jet (case LI) strong 
oscillations in the dilatation field developed near the centerline in the form of two-delta 
waves. They are evident in Fig. 2a in the region near the polar axis starting with the 
streamwise (x~15R 0 ) location where ring vortices begin to be shed from the jet. These point- 
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to point oscillations in dilatation are of the order of 1.5R(/U, which is more than one order of 
magnitude higher that the values of 0.1 R/U associated with the maximum dilatation inside 
the coherent structures of the jet. These problems were also apparent in the vorticity fields 
near the centerline. For case LI the simulation did not diverge, but obviously the quality of 
the solution, especially the sound information that can be collected from these fields is poor. 
When used to calculate a turbulent jet we were not able to get a well-behaved solution on the 
relatively coarse mesh used in these simulations without substantial filtering of the solution to 
remove the two-delta waves (every couple of time steps), and the quality of the solution was 
poor. This shows that the robustness of this polar axis treatment deteriorates rapidly on 
coarse meshes and for high Reynolds numbers, when strong non-linear interactions are 
present in the flow. Several features of this method are suspected to be at the origin of the 
above mentioned problems. First, the derivatives in the directions corresponding to the 
Cartesian coordinates are evaluated using six-order explicit central-differences instead of 6 th 
order Pade schemes, while in the evaluation of the radial derivatives with compact differences 
numerical boundary conditions using one-sided differences of lower order have to be 
formulated. Another source of errors may be associated with the bias introduced in the 

solution by the arbitrary choice (Ng/2 possibilities) of the two perpendicular directions after 

which the flow equations are solved at the centerline. 

Both of these approximations are automatically removed when a staggered mesh in the 
radial direction is employed in simulation L2. However, our results using the technique 
proposed in [1] showed that in the absence of filtering the solution would diverge after 
approximately 1000 time steps. That was due to a very local instability situated at the first 
point off the polar axis at a certain streamwise location. The dilatation field remained fairly 
smooth at all other streamwise locations close to the polar axis in contrast to the results seen 
in case LI. However, when a 6 th order accurate filter was applied every 300 time steps to 
smooth the velocity field, the dilatation contours shown in Fig. 2b remained smooth and no 
unphysical oscillations were observed near the polar axis. The behavior of the solution in 
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this first test case seems to suggest that this method is not very robust for strong non-linear 
problems, which would limit their use for turbulent calculations, especially on fairly coarse 
grids. 

Finally, in case L3 where the new treatment using asymptotic series expansion was used 
at the polar axis, the dilatation fields are smooth (Fig. 3c) and the flow features are very 
similar to the results obtained with the method of Mohseni and Colonius, but no filtering was 
necessary to avoid numerical instabilities to form near the polar axis. 

Based on the results from the first test case, we focus our attention for turbulent 
calculations on the method proposed by Mohseni and Colonius [1] and the method based on 
series expansions. Filtering the solution every 30 time steps despite being sufficient in 
keeping the solution from diverging, does not remove the unphysical spurious oscillations in 
the dilatation field that are observed very close to the polar axis at all streamwise locations for 
x>15R 0 (Fig. 3a). In fact the maximum values of the dilatation in these elongated structures 
very close to the polar axis are around 1.2R/U, which is significantly higher than the levels 
of dilatation recorded in the rest of the computational domain. That is true even if the filter is 
applied more often. One should point out that this amount of filtering was found to have 
non-negligible effects if sound sources were calculated from these fields. Thus, even though 
the method proposed in [1] gives comparable results with the series-expansions based 

technique for laminar forced jets, in turbulent simulations, the robustness and accuracy of the 

4 

method in [1] appears to be inferior. The turbulent case is much more difficult to be handled 
because the level of non-linear flow interactions is much higher and jet structures are passing 
through the centerline. Convergence of the solution on relatively coarse meshes can be 
obtained in our simulations only by filtering the solution when the method of Mohseni and 
Colonius [1] is used. 

In simulation T2, patches of high dilatation are observed in a region of about lOR^ in the 
streamwise direction starting at the end of the inviscid core. The position of these patches is 
not very close to the polar axis (Fig. 3b), except times when they are convected by the mean 
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jet motion through the polar axis. The maximum absolute value of the dilatation field in the 
turbulent region is at all times around O^R/U, including at the polar axis. The short waves 

that are seen at the top of Fig. 3b are rather related to the high aspect ratio Ax/Ar of the 

present grid, and are not a consequence of centerline instabilities. Away from the centerline 
the dilatation field appears smoother in case T1 compared to T2, but this is just an effect due 
to the use of the filtering in Tl. Even so, spurious waves or unrealistically high values of the 
dilatation (Fig. 3a - case Tl) or voiticity at the first 2-3 points off the centerline are not 
observed. 

The level of the unphysical spurious dilatation oscillations present in the domain close to 
the polar axis is greatly reduced compared to the case in which the method of Mohseni and 
Colonius [1] was used, while filtering was practically eliminated. The 6 th order filter was 
applied every 500 time steps to eliminate the high-frequency waves that form because of the 
high grid-stretching ratio. Overall there is a significant increase in the robustness for the 
method using asymptotic series expansions compared to the other two approaches. 
Moreover, because we are primarily interested in extracting sound information from the near 
fields, elimination of filtering that may compromise substantially the quality of the sound data 
is a key factor in assessing the different methods. 


4. SUMMARY 

In this paper a general method for handling the singularities arising at the polar axis of 
the governing flow equations in cylindrical coordinates based on power series expansions in 
the radial direction and Fourier series in the azimuthal direction was presented. Using the 
most general form of these series expansions a new set of equations at the polar axis was 
derived by calculating the limit of the various operators appearing in the governing system of 
equations (Eq. 1). The new method is computationally easy to implement and is less 
expensive than solving the equations in Cartesian coordinates at the centerline. The method 
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was proven to be well-defined in the sense that for any single-valued variable at the polar axis 
(scalar quantity or streamwise velocity) in the RHS of the Navier-Stokes equations only the 
m=0 mode was shown to have a non-zero coefficient in the limit of r— >0, despite the fact that 
individual terms had non-zero coefficients for modes ±2 and ±3. The momentum equations 

for u r and u 0 are coupled in the limit r— >0, and only modes m=± 1 are left in the RHS of the 

two momentum equations. Again this is consistent with the asymptotic behavior of the radial 
and azimuthal velocities at the origin (u r (0 + n/2) = u 0 (9)). 

The present algorithm avoids the loss in accuracy at the centerline where most of the 
finite-volume and finite-differences methods use some kind of one-sided differences to 
approximate the operators at the centerline. This is because in the context of the present 
method we were able to calculate the radial derivatives at the origin with the same order of 
accuracy as in the rest of the domain (6 th order compact differences) by using a domain 
mapping such as the points on the polar axis become regular interior points. As this is the 
only information needed to calculate the coefficients of the newly derived polar equations, 
there is no loss in the overall accuracy of the method. The coefficients in the asymptotic 
expansions are unique so we do not have to choose any arbitrary directions as is the case 
when the equations are solved in Cartesian coordinates at origin. One of the advantages is 
that the first point off the axis where the various terms with singular behavior have to be 

evaluated will be situated at Ar as opposed to Ar/2 (as is the case in the method of Mohseni 

and Colonius [1]), which may avoid an important source of errors or instabilities (e.g., due 
to the nonlinear interactions of eddies near the polar axis in a DNS or LES simulation). 

The robustness of the proposed approach was tested successfully by comparing results 
for a deterministically forced jet with similar calculations using same numerical method but 
with the equations solved at the centerline in Cartesian coordinates, or without any points at 
r=0 using the method described in [1], Finally, the present technique was shown to give 
improved results compared to the method of Mohseni and Colonius [1] for the simulation of 
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a compressible jet using LES on relatively coarse meshes, in which the flow interactions 
taking place near the polar axis are very important, and where a flow solution that is not 
contaminated by numerical artifacts originating at the singularity axis is required to accurately 
calculate the jet noise sources. 

One should point out that the algorithm presented here to deal with the polar-axis 
singularities is not restricted to the use of the present numerical method, that uses sixth-order 
Pade schemes to evaluate the derivatives in the radial and streamwise directions, or to solving 
the compressible flow equations. Rather it can be adapted in a straightforward way to handle 
the system of equations governing magneto-hydrodynamics, combustion, etc. As spherical 
coordinates are locally cylindrical coordinates near the two singularity axes, the above 
method is applicable to numerical methods that solve the flow equations in spherical 
coordinates. In this case the governing equations are different but the derivation of the set of 
equations valid at the two singularity axes is similar. 
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APPENDIX 

The mass conservation, momentum and equations for the compressible flow equations 
in cylindrical coordinates are: 


dp _ J dpu x + 1 drpUf 

dt v dx r dr 


+ ldpue 
r d0 


(Al) 


dpu x _ ( dpu x u x 1 drpu x u r 1 dptyV ) dp 

dt V dx r dr r d0 J dx x 


(A2) 
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The irreversible viscous dissipation, 'F, is given by: 
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The viscous terms appearing in the momentum equations are: 


V — !?Jxx | 1 | 1 9t x9 

* dx r dr r 90 

V — ^xr | ^ | j_ Tq9 

r dx r 9r r 90 r 

V — t 1 drigr 1 9r ee 
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while the viscous flux term in the energy equation is: 


9xvPr 9xJ r9r[ PrOrJ rOOVrPrOoJ 


(A8) 
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In Eq. (A8) T is the temperature, Pr is the Prandtl number and p = p(T) is the shear 
viscosity. The viscous stresses in cylindrical coordinates are: 



The coefficient A in (A9) is defined as: 


A = Mb_2 

p 3 


(A10) 


where p B is the bulk viscosity (p B /p = 0.6 for air according to [21]. The dilatation of the 
velocity field, 0 , is given by: 


0 = 


du 2L+ au L 

dx dr 



u r + 


due ) 

dej 


(All) 
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In the following the series expansions coefficients introduced in Eq. (3) for a scalar variable s 
are a-f.Pjftfor u x , s=x, for pu x , s=px), while those introduced in Eq. (4) for multi-valued 

variables (u r ,u 0 ) are A^.B^ with s=r and s=0. The derivative of the viscosity p(T) is 
calculated as: 


dp. _ dp dT 
dxj dT dxj 


= n( T) 


dT 

dx ; 


(A14) 


If the appropriate series expansions are introduced in Eqs. (Al) through (A8) and (A1 1) and 
the limit r-^0 is taken, a new set of equations valid at the polar axis is obtained. 
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In deducing the final form of Eqs. (Bl) through (B8) repeated use was made of Eq. (11). 
Eqs. (B3) and (B7) were obtained by using the expressions for the radial momentum 
equations and its viscous term (the first equation in (A3) and (A7)). As we pointed out in the 
paper the expressions obtained using the azimuthal momentum equation and its 
corresponding viscous term (the second equation in (A3) and (A7)) are redundant in the 
sense that the two independent equations corresponding to the m=l and m=-l modes will be 
exactly the same as the ones deduced using only the expressions of the radial momentum 
equation and its viscous term. Nevertheless, if this method will be applied to another system 
of equations, doing the calculations for both the radial and azimuthal directions may serve as 
a check of the calculation. For the LES simulations there are additional terms arising from 
the sub-grid scale modeling. For instance, the new term in the energy equation is: 
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where f> = pv t /Pr t , v, is the eddy viscosity and Pr, is the turbulent Prandtl number. As 6 

is a function of space that cannot be related directly to the temperature as was p, series 

expansions should be used for this coefficient when the limit is taken in Eq. (A 15). The final 
form will be: 
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dx 
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Note that 6 is evaluated directly in physical space using Lilly’s contraction in the present 
implementation of the dynamical model (as opposed of calculating first Pr,), and its series 

expansion is calculated in the same way as for the eddy viscosity v,. 
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FIGURES 


Fig.1 Total vorticity contours. Vorticity shown using 18 equally spaced contours between 
0 and 211/Rfl. (a) Laminar forced jet; (b) Turbulent jet 

Fig. 2 Dilatation contours in the forced jet. (a) Cartesian Coordinates; (b) Method of 
Mohseni and Colonius; (c) Series expansions treatment . 

Fig. 3 Dilatation contours in an area situated after the end of the potential core for the 
turbulent jet. Dilatation fields are shown with 18 equally spaced contours between 
-0.08U/R 0 and 0.08U/R 0 . (a) Method of Mohseni and Colonius; (b) Series expansions 


treatment 










